Multistability and Self-Organization in Disordered SQUID Metamaterials 
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Planar arrays of magnetoinductively coupled rf SQUIDs (Superconducting Quantum Interference 
Devices) belong to the emergent class of superconducting metamaterials that encompass the Joseph- 
son effect. These SQUID-based metamaterials acquire their electromagnetic properties from the 
resonant characteristics of their constitutive elements, i.e., the individual rf SQUIDs. In its simplest 
version, an rf SQUID consists of a superconducting ring interrupted by a Josephson junction. We 
investigate the response of a two-dimensional rf SQUID metamaterial to frequency variation of an 
externally applied alternating magnetic field in the presence of disorder arising from critical current 
fluctuations of the Josephson elements; in effect, the resonance frequencies of individual SQUIDs are 
distributed randomly around a mean value. Bistability is observed in the total current-frequency 
curves both in ordered and disordered SQUID metamaterials; moreover, bistability is favoured by 
disorder through the improvement of synchronization between SQUID oscillators. Relatively weak 
disorder widens significantly the bistability region by helping the system to self-organize itself and 
leads to nearly homogeneous states that change smoothly with varying driving frequency. Also, the 
total current of the metamaterial is enhanced compared with that of uncoupled SQUIDs, through 
the synergetic action of coupling and synchronization. The existence of simultaneously stable states 
that provide either high or low total current, allows the metamaterial to exhibit different magnetic 
responses that correspond to different values of the magnetic permeability, and provide either high 
or low total current, allows the metamaterial to exhibit two different magnetic responses that corre- 
spond to different values of its magnetic permeability. At low power of the incident field, high-current 
states exhibit extreme diamagnetic properties corresponding to negative magnetic permeability in 
a narrow frequency region. 



I. INTRODUCTION 

Advances in theory and nanofabrication techniques 
have opened many opportunities for researchers to cre- 
ate artificialy structured, composite media that exhibit 
extraordinary properties. The metamaterials (MMs) are 
perhaps the most representative class of materials of 
this type, which, among other fascinating properties, ex- 
hibit negative refractive index and optical magnetism 1 - - — . 
High-frequency magnetism, in particular, exhibited by 
the magnetic metamaterials, is considered one of the 'for- 
bidden fruits' in the Tree of Knowledge that has been 
brought forth by metamaterial research!. The unique 
properties of MMs are particularly well suited for novel 
devices like hyperlenses, which surpass the diffraction 
limits, and optical cloaks of invisibility 6 . Furthermore, 
they can form a material base for other functional devices 
with tuning and switching capabilities 4,7 . The key ele- 
ment for the construction of MMs has customarily been 
the split- ring resonator (SRR), a subwavelength "par- 
ticle" which is effectivelly a kind of an artificial "mag- 
netic atom"—. In its simplest version it is just a highly 
conducting ring with a slit that can be regarded as an 
inductive-capacitive resonant oscillator. SRRs become 
nonlinear and therefore tunable with the insertion of 
an electronic component (e.g., a diode) in their slits- 
However, metallic SRRs suffer from high ohmic losses 
that place a strict limit on the performance of SRR- 



based metamaterials, either in the linear or the nonlinear 
regime, and hamper their use in novel devices. The in- 
corporation of active constituents in metamaterials that 
provide gain through external energy sourses has been 
recognized as a promising technique for compensating 
losses^. On the other hand, the replacement of the met- 
alic elements with superconducting ones, provides both 
loss reduction and wideband tuneabilityii; the latter be- 
cause of the extreme sensitivity of the superconducting 
state to external stimuli 7,11 . Tunability of superconduct- 
ing metamaterial properties by varying the temperature 
or an externally applied magnetic field have been recently 
demonstrated^—. 

Going a step beyond, the metalic metamaterial ele- 
ments can be replaced by rf SQUIDs (rf Superconducting 
QUantum Interference Devices), creating thus SQUID- 
based metamaterial s 17 ' 18 . The rf SQUID, as shown in 
figure 1(a), consists of a superconducting ring interrupted 
by a Josephson junction (JJ) and demonstrates both re- 
duced losses and strong nonlinearities due to the Joseph- 
son elemen t 19 ' 20 . It constitutes the direct superconduct- 
ing analogue of a nonlinear SRR, that plays the role of 
the 'magnetic atom' in superconducting metamaterials in 
a way similar to that of the SRR for conventional (met- 
alic) metamaterial a 17 ' 21 ' 22 . However, currents and volt- 
ages in SQUID elements are determined by the celebrated 
Josephson relations 2 ^. The feasibility of using supercon- 
ducting circuits with Josephson junctions as basic ele- 



merits for the construction of superconducting thin-film 
metamaterials has been recently demonstrated 2 ^. Non- 
linearity and discreteness in SQUID-based metamaterials 
may also lead in the generation of nonlinear excitations 
in the form of discrete breather s 18 ' 21 ' 22 , time-periodic 
and spatially localized modes that change locally their 
magnetic response. Numerous types of SQUIDs have 
been investigated since its discovery, while it has found 
several technological application s 25 ' 26 . Recent advances 
that led to nano-SQUIDs makes possible the fabrica- 
tion of SQUID metamaterials at the nanoscale 2 ^. The 
use of SQUID arrays in dc current sensors 2 ^., filter o 29 ' 30 , 
magnetometers 3 "!,, amplifiers^—, radiation detectors 3 - 5 -, 
flux-to- voltage converters^ 6 -, as well as in rapid single flux 
quantum (RFSQ) electronics 3 -^, has been suggested and 
realized in the past. However, in most of these works 
the SQUIDs in the arrays were actually directly coupled 
through conducting paths. The SQUID-based metama- 
terial suggested in reference 1 ^ and investigated further 
in the present work relies on the magnetic coupling of its 
elements through dipole-dipole forces due to the mutual 
inductance between SQUIDs. 
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FIG. 1: (color online) (a) Schematic drawing of an rf SQUID 
in a perpendicular time-dependent magnetic field H(t). (b) 
Equivalent electrical circuit (RCSJ model) for a single rf 
SQUID driven by a flux source & ex t- 

Moreover, at low (sub-Kelvin) temperatures provide 
access to the quantum regime, where rf SQUIDs can 
be manipulated as flux and phase qubit o 38 ' 39 , the ba- 
sis element for quantum computation. Inductively cou- 
pled SQUID flux qubits can be used for the realiza- 
tion of quantum gates^S, while larger arrays of induc- 
tively coupled SQUID flux qubits have been proposed as 
scalable systems for adiabatic quantum computin g 41 ' 42 . 
Notably, amplification and squeezing of quantum noise 
has been recently achieved with a tunable SQUID-based 
metamaterial^ 3 -. 

In the present work we investigate the response of 
a two-dimensional (2D) rf SQUID metamaterial in the 
planar geometry with respect to frequency variation of 
an alternating magnetic field, focusing on the effect 
of quenched disorder through the SQUID parameter /3 
that determines the resonance frequency of individual 
SQUIDs. We are particularly interested in frequencies 
near resonance where multiple responses may appear; 
this issue is related to the existence or not of a bistability 



region in the current-frequency curve of the SQUID meta- 
material. Having calculated the response of the metama- 
terial to a given alternating filed, the effective magnetic 
permeability fi r can be determined by temporal and spa- 
tial averaging. Different responses, corresponding to dif- 
ferent simultaneously stable metamaterial states, result 
in different fj, r in the bistability region. In the next sec- 
tion we give a brief overview of the equations for a single 
SQUID, we introduce the equations for the flux dynam- 
ics in a 2D SQUID metamaterial and calculate its linear 
modes. In Section 3 we present numerical results for the 
maximum total current (devided by the total number of 
SQUIDs) as a function of the driving frequency. Such 
current-frequency curves are obtained both for ordered 
and disordered SQUID metamaterials, focusing primar- 
ily on the possibility of bistability. In Section 4 we discuss 
the effect of synchronization and present calculations of 
the relative magnetic permeability /i r . Section 5 contains 
the conclusions. 



II. DYNAMIC EQUATIONS AND LINEAR 
MODES 



In an ideal J J, the current-phase relation is of the form 
I = I c sm((f)j), where I c is the critical current of the J J 
and (j)j the Josephson phase. When driven by an ex- 
ternal magnetic field H(t), the induced (super)currents 
around the SQUID ring are determined by the celebrated 
Josephson relations^ 3 -. In the equivalent circuit picture, 
the resistively and capacitively shunted junction (RCSJ) 
model is frequently adopted to describe a real JJ. Thus, 
the equivalent lumped circuit for the rf SQUID in a mag- 
netic field with appropriate polarization comprises a flux 
source $> ex t in series with an inductance L and an ideal 
JJ, while the latter is shunted by a capacitor C and a 
resistor R [figure 1(b)]. Then, the dynamic equation for 
the flux threading the SQUID ring can be obtained by 
application of Kirkhhoff laws, as 
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= 0, (1) 



where <fr ex t is the external flux, $o is the magnetic flux 
quantum, and t is the temporal variable. The flux $ 
threading the SQUID ring is related to the Josephson 
phase through the flux quantization condition 



$ 



27m, 



(2) 



where n can be any integer. 

The rf SQUID is a highly nonlinear oscillator whose 
amplitude-frequency curves exhibit several peculiar fea- 
tures not seen in conventional inductive-capacitive os- 
cillators. For example, in a fine scale one may observe 
internal structure of increasing complexity that increases 
with increasing /S 4 ^. That fine structure of the amplitude- 
frequency curves can be reproduced numerically by inte- 
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FIG. 2: (color online) Maximum current amplitude i ma ,x as 
a function of the driving frequency fl for an rf SQUID with 
a = 0.002, tf> dc = 0, and (a) /3 = 1.27, cj> ac = 0.1; (b) /3 = 0.15, 
4>ac — 0.02. An enlargement of the bistability region is shown 
in the insets for each case. 



grating the dynamic equation (J) ) 18 i 22 . The SQUID res- 
onance can be tuned either by varying the amplitude of 
the alternating driving field or by varying the magni- 
tude of a static (dc) field threading the SQUID ring that 
creates a flux bias. The resonance shift due to nonlinear- 
ity has been actually observed in a Josephson paramet- 
ric amplifier driven by fields of different power levels 34 , 
while the shift with applied DC flux has been seen in 
high— T c rf SQUIDs 4 ^ and very recently in a low— T c rf 
SQUID in the linear regime 2 ^. Systematic measurements 
on microwave resonators comprising SQUID arrays are 
presented in reference d 34 ' 46 . For very low amplitude of 
the driving field (linear regime), the rf SQUID exhibits 
a resonant magnetic response at a particular frequency 
where luq 



ujsq = cjoVI + Pl, where luq — l/\/LC is the inductive- 
capacitive SQUID frequency and Pl is the SQUID pa- 
rameter 



p L = 2tt/3 = 2tt^. 



(3) 



The dynamic behavior of the rf SQUID has been stud- 
ied extensively for more than two decades both in the 
hysteretic (Pl > 1) and the non-hysteretic regimes, usu- 
ally under an external flux field of the form 



$ext = *d 



$ ac cos(wi), 



(4) 



where u is the driving frequency. The first and second 
term on the right-hand-side of the earlier equation corre- 
spond to the fluxes due to the presence of a constant (dc) 
and an alternating (ac) spatially uniform magnetic field, 



respectively. Typical current amplitude-frequency curves 
for a single SQUID are shown in figure 2 (for <&dc = 0). 
Equation ([TJ is formally equivalent to that of a massive 
particle in a tilted washboard potential (figure 3) 



U SQ - g 



$ P 



2L 



Ej cos 



2lr 7T 
$o 



(5) 



with Ej = I c $o/(2tt) being the Josephson energy. While 
for Pl < 1 there the potential has a single minimum, 
it aquires more and more local minima as Pl increases 
above unity. Moreover, applied dc flux moves the loca- 
tion of both the local and the global minima. 




FIG. 3: (color online) Potential curves as a function of the 
flux threading the SQUID ring, (a) For a non-hysteretic 
SQUID with /3 L ~ 0.75 < 1 and <j> dc = $ dc /$ = 
(black-solid curve); 0.5 (red-dashed curve); 1.0 (green-dotted 
curve), (b) For a hysteretic SQUID with /3l — 8 > 1 and 
(j>dc = &dc/&o — (black-solid curve); 0.5 (red-dashed curve); 
1.0 (green-dotted curve), (c) For $d c = and Pl — 0.5 < 1 
(black-solid curve); 1.5 (red-dashed curve); 2.5 (green-dotted 



Consider a planar array comprising identical rf 
SQUIDs arranged in a tetragonal N x x N y lattice, that is 
placed in a spatially uniform, alternating magnetic field 
directed perpendicularly to the plane of the SQUIDs. 
The flux threading each SQUID ring induces a current 
with both normal and superconducting components that 
generates its own magnetic field. The induced fields cou- 
ple the SQUIDs to each other through magnetic dipole- 
dipole interactions. The strength of this magnetoinduc- 
tive coupling falls off approximatelly as the inverse-cube 
of the distance between SQUIDs. If the distance between 
neighboring SQUIDs is such that they are weakly cou- 
pled, then next-nearest and more distant neighbor cou- 
pling can be neglected. In that case, we need to take 
into account nearest-neighbor interaction only between 



SQUIDs, and the flux threading (n,m)-th SQUID of 
the array is given by 

®n,m — ^ext + L [I n ,m + A x (v n _i iTO + I n +l,m) 

+ ^y(In,m-l + In,m+\)\ , (6) 

where n — 1,...,N X , m = l,...,N y , I n , m is the total 
current induced in the (n,m)— th SQUID, and X x ^ y = 
M x ^y/L are the magnetic coupling constants between 
neighboring SQUIDs in the x and y directions, respec- 
tively. The values of the M x and M y are negative since 
the magnetic field generated by the induced current in a 
SQUID crosses the neighboring SQUID in the opposite 
direction. By adopting the resistively and capacitively 
shunted junction (RCSJ) model, the current I n , m is 



d 2 $„, m 1 d$> n ,m , T {<, $n,m 

= C -TT 1 r J Mr sin 2tt 



dt 2 



R dt 
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(7) 



Then, following the procedure of reference^ and neglect- 
ing terms of order A^A^, A 2 , \ x , etc., we get 
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-A x ($n.-l,m + $n+l,m) — A tf ($n,m-1 + 3?rt,m+l) 
= [1 - 2(A X + \y)]$ext- 



(8) 



In the absence of losses (7 = 0), the earlier equations can 
be obtained from the Hamiltonian function 
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and 



C- 



d*,. 



df 



is the canonical variable conjugate to $ n . m , and repre- 
sents the charge accumulating across the capacitance of 
the JJ of each rf SQUID. The Hamiltonian function ([9]) is 
the weak coupling version of that proposed in the context 
of quantum computation 41 . Using the relations 
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and equations ([3]), (Tl"4"l) . equations © are normalized to 

<^n,m + 7<k,m + 4>n,m + P sin(27T0„ im ) 
— A x (0n-l,m + 0nH-l,m) — \{<fin,m-l + 0n,m+l) 

= &//> (12) 



3 

2 
1 

-l 

-2 
-3 



FIG. 4: (color online) Contours of the linear dispersion Q^ on 
the k x — K y plane for a two-dimensional rf SQUID array, with 
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-0.014 and /3 = 0.15. 



where the overdots denote differentiation with respect to 
the normalized time t, 



fieff = [1 _ 2 (^£ + Xy)]<t>ext, 

is the effective driving field, and 
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(13) 



(14) 



is the loss coefficient of individual SQUIDs, that actually 
represents all of the dissipation coupled to each rf SQUID 
and may also include radiative losses^. 

SQUID metamaterials support magnetoinductive flux- 
waves^, just like conventional metamaterials comprising 
metallic elements (i.e., split-ring resonators) 48 . The fre- 
quency dispersion for small amplitude flux waves is ob- 
tained by the substitution of <f> = A exp[i(K x n + K, y m — 
fir)], into the linearized equation dill) without losses and 
external field (7 = 0, 4> ex t = 0) 



n 



Jl+ f3 L ~ 2(\ x COS K X + \y COS Ky), (15) 

(1"J where ft = ui/coo and K x>y — d x ^ y k x . y are the normal- 



ized wavevector components, with k x (k y ) and d x (d y ) 
being the wavevector component and center-to-center 
distance between neighboring SQUIDs in x— direction 
( y— direction), respectively. 



III. CURRENT - FREQUENCY CURVES 

In the following, equations \\2\ are implemented with 
the boundary condition (except otherwise stated) 

<h,m(r) = 0JV a ,+l,m(l") = 4>n,o(r) = 4>n,N v + l(r) = 0,(16) 

for n = 1 , . . . , N x , m = 1, . .., N y that account for the ter- 
mination of the structure in a finite system. The set 



of dynamic equations (|16[) are integrated in time with a 
standard 4th order Runge-Kutta algorithm. For a sin- 
gle SQUID, the current denoted by i m ax is just the am- 
plitude of the induced current. In the case of ordered 
SQUID arrays we use the same notation for the maximum 

total current, i.e., i max = max Ijjjj^ Y, n ,m in,m(t)\, 
where n = 1 , . . . , N x and m — 1 , . . . , N y . For disordered 
arrays, the brackets < ... > indicate averaging of the 
maximum total current of the array over the number of 
differnt realizations, n r . In order to trace the current- 
frequency curves for ordered (i ma x — ^) or disordered 
(< i max > —Q) arrays, we start the system with zero ini- 
tial conditions and start integrating with a low (high) fre- 
quency until a steady state is reached. Subsequently the 
frequency is increased (decreased) by a small amount and 
the equations are again integrated until a steady state is 
reached, and so on. In each frequency step (except the 
first one, where we use zeros) we use as initial condition 
the steady state solution obtained in the previous step. 
The bistability properties of single SQUID oscillators (see 
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FIG. 5: (color online) Maximum total current amplitude imax 
as a function of the driving frequency Q for two-dimensional 
N x x N y rf SQUID arrays with a = 0.002, /3 = 1.27, <f> dc = 0, 
4>ac — 0.1, N x = N y — 20, and (a) periodic boundary condi- 
tions; (b) free-end boundary conditions, starting with differ- 
ent initializations. The black dotted lines indicate the corre- 
sponding imax vs. Q curves for a single rf SQUID. The maxi- 
mum total current of the arrays has been divided by the total 
number of rf SQUIDs N x x N y to facilitate the comparison. 



figure 2) are also seen in the arrays as well. We assume 
a moderate size array with N x — N y = 20, for which 
the coupling between SQUIDs is isotropic and weak, i.e., 
A x = X y = —0.014. Typical current-frequency curves are 
shown in figure 5, where the maximum of the total cur- 
rent, divided by the total number of SQUIDs in the array, 



is displayed as function of the frequency Q of an alternat- 
ing flux (dc flux is set to zero). In this figure, the value of 
(3— parameter has been selected so that hysteretic effects 
are rather strong (/3l — 8). We observe that bistability 
appears in a frequency region of significant width. The 
corresponding curves for a single SQUID are also shown 
for comparison. In figure 5(a), where periodic boundary 
conditions have been employed, we observe that although 
the bistability region for the array is narrower than that 
for a single SQUID, the maximum current per SQUID 
is slightly larger than that for a single SQUID. In the 
case of periodic boundary conditions, the size of the ar- 
ray does not affect those results; current-frequency curves 
for larger arrays with N x = N y = 40 and N x — N y = 80 
(not shown) are practically identical to these shown in 
figure 5(a). In figure 5(b) and the rest of the paper free- 
end boundary conditions [equations (|16[) ] have been used. 
In this case, the current-frequency curves are very sen- 
sitive to the initial conditions, the frequency step, and 
other parameters. The curves in figure 5(b) shown in 
different colours (red and green) correspond to different 
initializations of the same system. 
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FIG. 6: (color online) Maximum total current amplitude i max 
as a function of the driving frequency fi for two-dimensional 
N x x N y rf SQUID arrays with a = 0.002, = 0.15, cj> dc = 0, 
0ac = 0.02, and (a) N x = N y = 20; (b) N x = N y = 40. 
The black dotted lines indicate the corresponding i max vs. fi 
curves for a single rf SQUID. The maximum total current of 
the arrays has been divided by the total number of rf SQUIDs 
N x x N y to facilitate the comparison. Free-end boundary 
conditions that account for the termination of the structure 
have been used. 

The parts of current-frequency curves that are close to 
those for the single SQUID are formed by almost ho- 
mogeneous states, i.e., states with all SQUIDs in the 
high-current or low-current state. Homogeneous states 
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FIG. 7: (color online) Maximum current amplitude averaged 
over Nr = 30 realizations od disorder, < i ma x >, as a func- 
tion of the driving frequency £1 for a SQUID metamaterial 
with a = 0.002, = 1.27, (f> ac = 0.1, dc = 0, and (a) 
j3 = 1.27 ±0.01; (b) p = 1.27 ±0.05; (c) fi = 1.27 ±0.1. 
The corresponding curves for the same SQUID metamaterial 
without disorder are shown in black-dotted lines. 
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FIG. 8: (color online) Maximum current amplitude averaged 
over Nr = 30 realizations od disorder, < i ma x >, as a func- 
tion of the driving frequency fi for a SQUID metamaterial 
with a = 0.002, /3 = 0.15, 4> ac = 0.02, 4> dc = 0, and (a) 
/3 = 1.27 ±0.01; (b) /3 = 1.27 ±0.05; (c) /3 = 1.27 ±0.1. 
The corresponding curves for the same SQUID metamaterial 
without disorder are shown in black-dotted lines. 



are formed easier in the periodic systems [figure 5(a)]; 
only a small part survives in the case of free bound- 
ary conditions [figure 5(b)] that is close to the current- 
frequency curve of the single SQUID. In the latter figure 
we also observe the formation of small steps for which the 
corresponding solutions may be characterized as 'mixed 
states', that are formed by a certain number of SQUIDs 
in the high-current state while all the others are in the 
low-current state. In figure 6, the corresponding curves 
for SQUIDs with /3l — 1 (P = 0.15) are shown. A com- 
parison with the corresponding curves for a single SQUID 
(shown as dotted black lines) indicates that the frequency 
region where bistability appears are nearly of the same 
width. 

The assumption of SQUID-based metamaterials com- 
prising identical elements is certainly not realistic. 
Therefore, we consider disordered SQUID arrays in which 
the parameter /3 varies randomly within a particular 
range of values around a mean value. The SQUID pa- 
rameter /?, that depends on the critical current of the 
Josephson junctions, determines also the resonance fre- 
quency of individual SQUIDs. We have calculated the 
maximum current-frequency curves of the same SQUID 
metamaterials as before, taking into account the distribu- 
tion of the natural frequencies of individual SQUIDs. For 
obtaining reliable results, we have taken statistical aver- 
ages over many realizations, ur, of disorder. Remark- 
ably, the calculations reveal that weak disorder strongly 



favours bistability, as it is observed in figures 7 and 8 for 
mean values of j3 = 1.27 and 0.15, respectively. In these 
figures, as we go from (a) to (c) j3 fluctuates by ±0.01, 
±0.05, and ±0.1, so that the relative disorder strength 
is much larger in figure 8. In all cases shown, either 
exhibiting weak or strong disorder, the stability of the 
nearly homogeneous states with high current increases 
considerably with respect to that for the corresponding 
ordered SQUID metamaterials. This is actually reflected 
in the widening of the bistability regions of the current 
frequency curves. We should also note that the bistabil- 
ity region gradually shrinks with increasing strength of 
disorder. That shrinking occurs more rapidly for /3j, ~ 1 
(figure 8) since the relative /3 variation is larger in this 
case. 



IV. SYNCHRONIZATION AND 
MULTIRESPONSE 

In order to ensure that the maximum current- 
frequency curves presented in figures 7 and 8 corre- 
spond to homogeneous states, we define and calculate 
a (Kuramoto-type) synchronization parameter, as 



# = 



1 >p e 2ni4 

N a N„ ^ 

n.m 



(17) 



v y 



where the brackets denote averaging both in time (i.e., 
in one oscillation period) and the number of realizations 
of disorder tir. The absolute value of "J/ quantifies the 
degree of synchronization; |\&| may vary between and 
1, corresponding to completely asynchronous and syn- 
chronized states, respectively. The calculated values for 
(parts of) the maximum current-frequency curves in fig- 
ure 7 are shown in figure 9 for J7 varying in both direc- 
tions. The bistability regions are shown in this figure as 
green-dotted vertical lines, for reference. For very weak 
disorder [figure 9(a)], |\P| remains close to unity for most 
of the frequency interval shown. However, there is a nar- 
row region at low frequencies where synchronization, and 
therefore complete homogeneity breaks down, when the 
high current solution loses its stability. In that case all 
the SQUIDs change their state towards a lower maxi- 
mum current state. During the proccess, the phase dif- 
ferences of the SQUID fluxes and the driving field lock 
in randomly selected values, resulting in a low current 
and only partially synchronized state. However, the two 
low current states that can be distinguish in this fre- 
quency region, i.e., the synchronized one obtained with 
increasing frequency and the partially synchronized one 
obtained with decreasing frequency, provide almost the 
same maximum current. With increasing the strength 
of the disorder [figure 9(b)] the bistability region shinks 
while the same effect as in figure 9(a) is observed in a 
wider frequency interval. Moreover, the high current 
states are not completely synchronized in the bistabil- 
ity region, since \^\ is slightly less than unity. The ef- 
fect can be seen more clearly by further increasing the 
strength of the disorder as in figure 9(c), where |\&| is 
clearly less than unity in the bistability region indicating 
partial synchronization. 

The magnetic response of the SQUID metamaterial at 
a particular state can be calculated in terms of the mag- 
netization along the lines given in reference d 17 ' 24 . As- 
suming a tetragonal unit cell with d x — d y = d and a 
squared SQUID area of side a, the magnetization is 
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the spatially and temporally averaged current, and D is 
a length related to the cavity where the metamaterial 
is placed 24 . Using fundamental relations of electromag- 
netism we write the relative magnetic permeability as 
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FIG. 9: (color online) Fhe magnitude of the synchronization 
parameter |\&| as a function of the driving frequency Q in the 
bistability region, for a N x x N y = 20 x 20 SQUID metama- 
terial with a = 0.002, /3 = 1.27 ± 0.01, ac = 0.1, (f>dc = 0, 
Nr = 30 realizations of disorder, and (a) ft = 1.27 ± 0.01; 
(b) ft = 1.27 ± 0.05; ft = 1.27 ± 0.1. The arrows indicate the 
direction of frequency variation while the green dotted lines 
the corresponding bistability intervals. 
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FIG. 10: (color online) (a) Relative magnetic permeability 
fi r = m/Mo for the low and high maximum current states as 
a function of the driving frequency Q, for N x = N y — 20, 
a = 0.002, = 1.27, 4> a c = 0.1, and (j> dc = 0. Multiple-valued 
magnetic response is observed in the bistability region, (b) 
The corresponding maximum current-frequency curves. The 
arrows indicate the direction of frequency variation. 



where H is the intensity of a spatially uniform magnetic 
field applied perpendicularly to the SQUID metamaterial 
plane. The latter is related to the external flux to the 
SQUIDs as 



where (j,q is the magnetic permeability of the vacuum, 
and the brackets denote temporal averaging. Combining 
equations (fT8")) - PD|) , we get 
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FIG. 11: (color online) (a) Relative magnetic permeability 
fir — fi/fio for the low and high maximum current states as 
a function of the driving frequency Q, for N x = N y — 20, 
a = 0.002, p = 1.27, cj> a c = 0.001, and 4> dc = 0. Nega- 
tive fi r is observed in a narrow frequency band just above 
the resonance frequency, (b) The corresponding maximum 
current-frequency curves. 



1 

0.8 

o 

i 0.6 

0.4 
0.4 

8 
. £ 0.2 



1 1 


' 1 

1 


1 1 

(a)- 

l 


-- -N 


- \ 

i 

2 


1 

1 
3 


1 
(b)" 

1 
4 



Q 



FIG. 12: (color online) (a) Relative magnetic permeability 
fi r — fi/fMs for the low and high maximum current states as a 
function of the driving frequency Q,, for a disordered SQUID 
metamaterial, with N x = N y = 20, a = 0.002, fS = 1.27±0.1, 
4>ac — 0.1, and 4>dc = 0. Multiple-valued magnetic response 
is observed in the bistability region, (b) The corresponding 
maximum current-frequency curves. 



where k — tis §- £ - Jtq • For a rough estimation of the con- 
stant k assume that L ~ Hqol, where L is the SQUID 
inductance, and that D ~ d. Then, we have that 
k-/3(§) 3 . Using a = d/2 and /3 = 1.27 we get k ~ 0.16. 
We may then use the numerically calculated values of 
i n ,m and cf> ex t into equation (f2~Tj) to obtain fi r . Simultane- 
ously stable SQUID metamaterial states respond differ- 
ently to the external field and therefore exhibit different 
fi r . This can be seen clearly in figure 10, where the rel- 
ative permeability fi r has been calculated from equation 
(j2"Tj) . The SQUID metamaterial for the parameters used 
in figure 10 is diamagnetic for all frequencies; however, 
the diamagnetic response is stronger for the high current 
states in the bistability region. For a weaker driving field 
that provides a flux amplitude of 10 -3 the metamaterial 
is in the linear limit, as can be infered by inspection of 
the current-frequency curve shown in figure 11(b). The 
corresponding fi r as a function of the driving frequency 
f2 is again diamagnetic everywhere except close to the 
resonance, where strong variation of fx r occurs. For fre- 
quencies below (but very close to) the resonance at Q ~ 3 
the metamaterial becomes strongly paramagnetic. To the 
contrary for frequencies above (but very close to) the 
resonance the metamaterial becomes extremely diamag- 
netic, exhibiting negative fi r within a narrow frequency 
region. Note that in this case negative fi r would have 
also been obtained with a much smaller coefficient k. 

For a disordered SQUID metamaterial, that is rela- 
tively strongly driven (0 QC = 0.1), the relative permeabil- 
ity fi r changes only slightly. The most important effect 
observed in this case is the enlargement of the bistabil- 
ity interval (figure 12). For weakly driven, disordered 
SQUID metamaterial, increasing disorder results in de- 
creasing of the magnitude of fx r in the frequency region 
around resonance. This effect is illustrated in figures 13 
and 14, obtained for two different values of disorder; j3 
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FIG. 13: (color online) (a) Relative magnetic permeability 
fir — fl/flo for the low and high maximum current states as a 
function of the driving frequency Q,, for a disordered SQUID 
metamaterial, with N x = N y = 20, a = 0.002, [3 = 1.27±0.01, 
4>ac = 0.001, and <f>dc = 0. Negative fi r is observed in a narrow 
frequency band just above the resonance frequency, (b) The 
corresponding maximum current-frequency curves. 



fluctuates around its mean value by ±0.01 and ±0.1, re- 
spectively. In the frequency region around resonance, 
in particular, we observe in figure 13 that the dip cor- 
responding to negative fi r becomes shallower than that 
for the ordered SQUID metamaterial (figure 11). With 
further increasing disorder (figure 14), the negative fi r 
region disappears. For this particular set of parameters, 
the minimum of fi r touches the zero axis. 



V. CONCLUSIONS. 

We investigated numerically two-dimensional SQUID 
metamaterials driven by an alternating magnetic field. 




FIG. 14: (color online) (a) Relative magnetic permeability 
fi r — m/^o f° r the low and high maximum current states 
as a function of the driving frequency £2, for a disordered 
SQUID metamaterial, with N x = N y = 20, a = 0.002, /3 = 
1.27 ± 0.1, 4> ac = 0.001, and <f>dc — 0. Negative /^ r is not 
observed because of the relatively stong disorder, (b) The 
corresponding maximum current-frequency curves. 



We have calculated current-frequency curves both for or- 
dered and disordered metamaterials; in both cases we ob- 
served bistability regions created by almost homogeneous 
high and low current states. However, we observed that 
the presence of disorder widens significantly the bistabil- 
ity regions. The effect is rather strong for SQUID meta- 
materials comprising SQUIDs with either small or large 
Pl parameter. Remarkably, the homogeneity of high and 
low current states, that is quantified by the parameter "J/, 
persists also in the case of disorder up to a very high de- 
gree. These states deviate only slightly from practically 
complete homogeneity only in the case of strong disorder, 
as can be seen from calculations of the amplitude of ^>. 
This important result indicates that random variation of 
the SQUID parameters does not destroy bistability, that 
is crusial for applications that require bistable switch- 
ing properties, but instead stabilizes the system against 
modulational or other instabilities. 

This result is related to past work on disordered net- 
works of nonlinear oscillators where it was concluded that 
moderate disorder may enhance synchronization and sta- 
bilize the system against chao a 49 ' 50 . In the present con- 
text, synchronization of individual SQUIDs in the high or 
low current states results in high or low maximum total 
current for the metamaterial. This requires that (almost) 
all the SQUIDs are in phase. It could be natural to as- 



sume that the more nearly identical the elements, the 
better the synchronization will be. However, even in the 
ideal case of identical elements, the earlier assumption 
may not be true and the in phase state may be dynami- 
cally unstable. Then, synchronization is reduced and the 
SQUID metamaterial cannot remain in the high current 
state that is more sensitive to instability. This can be 
clearly observed in figures 5(b) and 6(b), where in most 
of the bistability region the metamaterial relaxes to par- 
tially synchronized states that provide significantly lower 
maximum total current. This type of disorder-assisted 
self-organization may also occur by introducing local dis- 
order in an array of otherwise identical oscillators, i.e., in 
the form of impuritie s 51 ' 52 . In this case, the impurities 
trigger a self-organizing process that brings the system 
to complete synchronization and suppression of chaotic 
behaviour. 

Having calculated numerically the response of the 
metamaterial to an alternating field with given frequency, 
we have calculated the magnetic permeability of the 
metamaterial for several illustrating cases both with and 
without disorder. While the expression for calculating 
the magnetic permeability is rather simple, there is some 
uncertainty about the value of the factor n. However, 
for a reasonable value of k we observe that the magnetic 
permeability can take negative values in a narrow fre- 
quency region above resonance for weakly driven SQUID 
metamaterials. In this case, increasing disorder results 
in weakening the negative response of the metamaterial; 
thus, for relatively strong disorder the response is not 
sufficient to make the magnetic permeability negative. 
For SQUID metamaterials exhibiting bistability, differ- 
ent magnetic permeabilities can be reached under the 
same conditions depending on their state. 
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